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ABSTRACT 

This  paper  presents  a  new  finite  difference  scheme  for  the  Stokes 
equations  and  incompressible  Navier-Stokes  equations  for  low  Reynold's 
number.  The  scheme  uses  the  primitive  variable  formulation  of  the  equations 
and  is  applicable  with  non-uniform  grids  and  non-rectangular  geometries. 
Several  other  methods  of  solving  the  Navier-Stokes  equations  are  also  examined 
in  this  paper  and  some  of  their  strengths  and  weaknesses  are  described. 
Computational  results  using  the  new  scheme  are  presented  for  the  Stokes 
equations  for  a  region  with  curved  boundaries  and  for  a  disc  with  polar 
coordinates.  The  results  show  the  method  to  be  second-order  accurate. 
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SIGNIFICANCE  AND  EXPLANATION 


The  incompressible  Navier-Stokes  equations  describe  the  flow  of  many 
common  fluids  such  as  water  or  air  at  low  speeds.  Thus  numerical  methods  for 
solving  these  equations  are  very  important  for  many  scientific  and  engineering 
applications.  In  this  paper  a  new  second-order  accurate  finite  difference 
method  for  solving  the  incompressible  Navier-Stokes  equations  is  presented. 

An  advantage  of  this  new  scheme  over  other  methods  is  that  it  can  be  applied 
with  non-rectangular  regions  or  non-uniform  grids. 

In  this  paper  several  other  methods  for  solving  the  Navier-Stokes 
equations  are  examined  and  some  of  their  strengths  and  weaknesses  are 
discussed. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


FINITE  DIFFERENCE  METHODS  FOR  THE 
STOKES  AND  NAVIER-STOKES  EQUATIONS 


John  C.  Strikwerda 


1.  INTRODUCTION.  In  this  paper  we  examine  several  common  methods  for 
solving  the  incompressible  Navier-Stokes  equations  by  finite  differences  and 
we  present  a  new  second-order  accurate  finite  difference  scheme  for  these 
equations.  This  new  scheme  is  designed  to  be  applied  with  non-uniform  grids 
and  non-orthogonal  coordinate  systems.  Numerical  experiments  with  the  Stokes 
equations  illustrate  the  versatility  and  accuracy  of  the  scheme. 

The  steady-state  Stokes  equations  on  a  domain  8  in  Rn  are  given  by 


(1.1) 


V2£  +  $p  -  f(x) 
$*u  *  g(x) 


and  the  steady-state  Navier-Stokes  equations  are 


(1.2) 


-R~VS  +  (u*$)u  +  $p  *  f ( x) 
$»u  *>  g(x) 


where  R  is  the  Reynolds  number.  He  will  consider  the  systems  (1.1)  and 

(1.2)  with  Dirichlet  boundary  conditions 

(1.3)  u(x)  *  b(x)  on  38  . 


A  necessary  condition  for  (1.1)  or  (1.2)  to  have  a  solution  is  that  the 
data  g(x)  and  b(x)  satisfy  the  integrability  condition 

(1.4)  /9s*/  b*n  , 

8  38 

where  n  is  the  outer  unit  normal  to  38.  For  the  mathematical  theory  of  the 
systems  (1.1)  and  (1.2)  we  refer  to  Ladyzhenskaya  (1963)  and  Teman  (1979). 

We  will  be  concerned  only  witlj  methods  that  solve  the  systems  (1.1)  and 
(1.2)  in  the  primitive  variables  u  and  p  and  not  with  methods  such  as  the 
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vorticity  and  stream-function  reformulation.  Also  our  methods  are  applicable 
in  two  or  three  dimensions  although  our  examples  will  be  only  in  two 
dimensions. 

We  emphasize  that  the  scheme  presented  here  is  designed  to  be  easily 
applicable  with  non- rectangular  geometries  and  non-uniform  grids.  The  vast 
majority  of  papers  on  the  numerical  solution  of  the  incompressible  Navier- 
Stokes  equations  limit  themselves  to  examples  using  rectangular  geometry  and 
uniform  grids.  By  way  of  contrast,  computations  with  the  compressible  Navier- 
Stokes  equations  routinely  use  non-rectangular  geometries  and  non-unifom 
grids. 

The  outline  of  the  remaining  sections  of  the  paper  is  as  follows.  In 
Section  2  we  discuss  the  strengths  and  weaknesses  of  some  common  approaches  to 
solving  the  systems  (1*1)  and  (1.2)  and  in  Section  3  we  discuss  finite 
difference  schemes  for  these  systems.  The  finite  difference  integrability 
condition  is  discussed  in  Section  4,  and  computational  results  are  given  in 
Section  5.  The  numerical  examples  of  Section  5  demonstrate  that  the  new 
scheme  given  here  can  be  used  to  give  second-order  accurate  solutions  to  the 
Stokes  equations  for  non-rectangular  geometries.  To  our  knowledge  no  other 
finite  difference  schemes  for  the  Stokes  or  incompressible  Navier-Stokes 
equations  in  the  primitive  variables  have  been  shorn  to  be  second-order 
accurate  for  non-rectangular  geometries.  Computations  using  the  new  scheme 
for  the  incompressible  Navier-Stokes  equations  are  currently  being  made  and 
will  be  reported  when  complete. 

2.  SOLUTION  TECHNIQUES.  In  this  section  we  review  some  approaches  to 
solving  the  Navier-Stokes  and  Stokes  equations  numerically.  Few  researchers 
have  treated  the  system  (1.2)  in  the  given  form,  most  have  altered  it  in  some 
way.  Before  examining  the  altered  forms  of  (1.2)  we  look  at  the  system  in  the 
given  form. 

The  Stokes  equations  (1.1)  and  the  Navier-Stokes  equations  (1.2)  are 
elliptic  systems  of  n  +  1  equations  in  n  +  1  dependent  variables.  The 
definition  of  an  elliptic  system,  as  given  by  Douglis  and  Nirenberg  (1957), 
requires  that  the  determinant  of  the  principle  symbol  of  the  system  not  vanish 
for  non-zero  values  of  dual  variables.  For  the  Navier-Stokes  equations  the 
determinant  of  the  principle  symbol  is 


which  is  non-zero  for  |€|  *  0.  Moreover,  since  the  determinant  is  a 
polynomial  of  degree  2n  in  the  variables  €  **  (€1,...,€  )  the  system 
requires  n  boundary  conditions  at  each  point  of  the  boundary  (Agmon,  Douglis 
and  Nirenberg  (1964)).  These  boundary  conditions  will  usually  be  Dirichlet  or 
Neumann  conditions  on  the  velocity  u. 


One 
(1.2)  is 


(2.2) 


of  the  most  common  ways  of  modifying  the  Navier-Stokes  equations 
to  replace  it  by  the  system 

-R  ^u  +  (u*$)u  +  $p  -  f  ( x ) 

V2p  -  $«f  +  R_Vg  -  l  u1  uj  -  u'tg 

i»  J  j  i 


The  last  equation  of  (2.2)  is  obtained  by  taking  the  divergence  of  the  first 
equation  of  (1.2)  and  then  using  the  last  equation  of  (1.2)  to  eliminate  the 
divergence  of  velocity.  The  system  (2.2)  has  the  advantage  over  (1.2)  in 
that,  when  discretized,  it  can  be  solved  using  standard  methods  for  inverting 
the  discrete  Laplacian.  However,  the  system  (2.2)  has  a  grave  disadvantage  in 
that  it  requires  n  +  1  boundary  conditions,  one  for  each  elliptic  equation, 
as  opposed  to  (1.2)  which  requires  only  n  boundary  conditions.  Thus  any 
attempt  to  solve  (1.2)  via  (2.2)  would  require  some  means  of  determining  the 
correct  additional  boundary  condition.  Without  the  correct  condition 
solutions  of  (2.2)  will  not  be  solutions  of  (1.2). 


Roache  (1972,  p.  194)  suggests  that  the  additional  boundary  condition  be 
given  by  the  normal  derivative  of  pressure  as  determined  by  the  first  equation 
of  (1.2)  or  (2.2)  evaluated  on  the  boundary.  This,  however,  is  not 
satisfactory  as  a  boundary  condition  since  it  is  not  independent  of  the  system 
of  differential  equations.  Roache's  suggestion  leaves  the  system  (2.2) 
underdetermined. 


Another  boundary  condition  which  is  commonly  used  along  boundaries 
corresponding  to  physical  surfaces  is  to  set  the  normal  derivative  of  the 
pressure  to  zero,  which  is  valid  in  the  limit  for  high  Reynolds  number  flow. 
With  this  boundary  condition  and  (1.3)  the  system  (2.2)  has  the  proper  number 
of  boundary  conditions,  however,  its  solutions  are  not  solutions  of  (1.2). 

As  one  would  expect,  the  methods  using  (2.2)  or  similar  systems  have 
difficulty  with  the  accuracy  of  the  pressure  field  and  with  satisfying  the 
incompressibilty  condition  on  the  velocities  (see  for  example  the  work  by 
Boney,  Hefner,  Hirsh,  and  Zoby  reported  in  Rubin  and  Harris  (1975)). 

The  above  mentioned  difficulties  are  seen  in  computations  with  the  time- 
dependent  Navier-Stokes  equations  as  well.  Roache  (1972)  has  a  discussion  of 
the  difficulties  of  obtaining  a  zero  divergence  for  the  velocity  field  when 
using  the  above  approach  for  time  dependent  flows  (see  also  Harlow  and  Welch 
(1964)). 

Because  of  these  difficulties,  it  seems  best  not  to  use  the  derived 
system  (2.2)  but  to  use  the  original  system  (1.2). 

Another  approach  to  solving  the  Navier-Stokes  equations  (1.2)  is  the 
artificial  compressibility  method.  The  basic  idea  of  this  method  is  to  solve 
a  time-dependent  system  of  equations,  whose  steady-state  solutions  solve 
(1.2),  until  a  steady  state  is  reached.  Methods  have  been  proposed  by  Chorin 
(1967)  and  Yanenko  (1967).  The  convergence  rate  of  these  methods  is  dependent 
on  the  choice  of  finite  difference  method  used  to  solve  the  system.  Moreover, 
as  will  be  discussed  in  Section  4,  it  may  happen  that  the  finite  difference 
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equations  do  not  have  a  steady-state  solution,  so  the  method  cannot 
converge.  Taylor  and  Ndefo  (1970)  reported  difficulty  in  getting  Yan. nko*s 
method  to  converge,  most  likely  because  there  was  no  solution. 

Another  common  method  is  to  use  the  "parabolized"  Navier-Stokes  equations 
in  which  the  second-derivatives  in  the  stream-wise  direction  are  removed. 
Because  of  its  limited  applicability  and  uncertain  justification  we  will  not 
discuss  this  method  here  except  to  note  that  often  an  analogue  of  (2.2)  is 
derived  and  thus  some  of  our  observations  on  (2.2)  also  apply  to  the 
parabolized  equaitons.  Raithby  and  Schneider  (1979)  discuss  these 
difficulties  for  three-dimensional  flow  problems. 

3.  FINITE  DIFFERENCE  SCHEMES.  In  this  section  we  discuss  the  staggered 
mesh  and  central  finite  difference  schemes  for  (1.1)  and  (1.2)  and  introduce  a 
new  scheme.  The  second-order  accurate  staggered  mesh  scheme  for  a  uniform 
cartesian  grid  assigns  the  values  of  each  of  the  velocity  components  and  the 
pressure  to  different  interlaced  grids.  In  two  dimensions  with  velocity 
components  u  and  v,  one  may  assign  values  of  u  to  grid  locations 

((i  +  -~)h, jh),  values  of  v  to  (ih,(j  +  ^)h),  and  values  of  p  to 

(ih,jh),  e.g.  Harlow  and  Welch  (1965),  Patankar  and  Spalding  (1972),  Raithby 
and  Schneider  (1979),  Brandt  and  Dinar  (1979).  This  method  works  very  well  as 
long  as  the  geometry  is  rectangular  and  the  grid  is  uniform.  Non-uniform 
grids  and  grid  mapping  techniques  cannot  be  conveniently  handled. 

The  staggered  mesh  schemes  also  have  some  difficulty  at  boundaries.  For 
example ,  when  both  velocity  components  are  specified  at  a  boundary  then  that 
velocity  component  whose  mesh  lines  do  not  lie  on  the  boundary  requires  some 
special  treatment. 

The  central  difference  scheme  on  a  uniform  rectangular  mesh  assigns 
values  of  all  the  variables  to  each  grid  point.  The  divergence  and  gradient 
operators  are  approximated  using  central  differences  and  the  Laplacian  is 
approximated  by  the  standard  five-point  discrete  Laplacian.  Central 
difference  schemes  have  been  used  by  Chorin  (1967,  1968)  in  time-dependent 
calculations. 

An  important  concept  for  finite  difference  schemes  for  elliptic  systems 
such  as  (1.1)  and  (1.2)  is  that  of  regularity  (see  Bube  and  Strikwerda  (1980), 
and  also  Frank  (1968),  Brandt  and  Dinar  (1979)).  Regular  schemes  give  rise  to 
regularity  estimates  analogous  to  those  in  the  theory  of  elliptic  systems  of 
differential  equations.  Solutions  to  regular  difference  schemes  will  in 
general  be  smoother  than  solutions  to  non-regular  schemes  and  also  will  be 
more  accurate  approximations  to  the  solutions  of  the  differential  equations. 

The  central  difference  scheme  is  non-regular  (Bube  and  Strikwerda 
(1980)),  which  results  in  non-smooth  solutions.  The  lack  of  smoothness  is 
most  noticeable  in  the  pressure.  The  staggered  mesh  scheme  is  regular.  The 
advantage  of  the  central  difference  scheme  is  that  it  is  easily  implemented 
with  non-uniform  grids  as  introduced  by  coordinate  changes. 

It  should  be  emphasized  that  none  of  the  difficulties  mentioned  above  are 
insurmountable.  Both  the  staggered  mesh  and  central  differencing  schemes  have 
been  used  and  often  quite  successfully.  However  we  will  consider  a  new  scheme 
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which  incorporates  both  regularity  and  ease  of  implementation  with  coordinate 
grid  mapping  techniques • 

Before  introducing  the  new  scheme  we  will  discuss  the  concept  of 
regularity  for  difference  schemes  as  given  in  Bube  and  Strikwerda  (1980).  A 
difference  operator  A  may  be  written  as 

Af (x)  -  l  a  (h,x)TUf(x)  , 

V  w 

where  TW  is  the  translation  operator  given  by 

TWf(xv)  -  ftx^) 

for  multi-indices  v  and  u* 

The  symbol  of  A  is  given  by 

a(h,x,C)  -  l  a  (h,x)eiW*C  . 

V  V 

For  example,  the  first-order  central  difference  operator  in  the  k-th 
coordinate  direction  has  symbol 

*k  ~«k 

5 - 5? - ihk  8in  S 

k 

and  the  standard  second-order  accurate  Laplacian  in  n-variables  has  the  symbol 


-  J,  i  s  • 

A  finite  difference  scheme  for  the  Stokes  equations  is  regular  elliptic  if  the 
determinant  of  the  matrix  of  symbols  of  the  scheme  vanishes  only  for  C  equal 
to  zero  modulo  2*.  For  the  Stokes  equations  with  central  differencing,  and 
Ax  ■  Ay  m  h,  this  determinant  is 


-2.  ,21.  ^  ,21, 

4h  (sin  -  C1  +  sin  —  ?2 


ih  'sin  C. 


.. -2,  ,21,  x  ,21 

4h  (sin  —  +  sin  —  C2) 

ih"'sin 


ih  'sin  C 


ih  'sin  C, 


-  4h"4(sin2  ~  C  +  sin2  ^  C  Hsta2^  +  sin2C2)  . 
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This  determinant  vanishes  for  the  dual  variables  ? 1  and  ?2  equal  to  *, 
and  thus  the  scheme  is  not  regular.  One  sees  that  the  non-regularity  comes 
from  the  form  of  differencing  used  for  the  gradient  and  divergence  terms.  Our 
new  scheme  is  a  modification  of  the  central  differencing  scheme  so  as  to  make 
the  scheme  regular. 


The  new  scheme  we  consider  will  be  called  the  regularized  central 
difference  scheme.  In  this  scheme  the  derivatives  of  pressure  are 
approximated  as 

“  6  p  -  ah?6  62  p 
kO  k  k-  k+* 

and  the  first  derivatives  of  the  velocity  in  the  divergence  equation  are 
approximated  as 


(3.3) 


JH*.  6 

3xk  k0 


uk  -  ch" 


2.  A2  * 

\A-U 


where  a  is  a  non-zero  constant  and  6  ,  6  and  6  are  the  centered, 

ItO  v— 

forward,  and  backward  divided  differences,  respectively.  The  Laplacian  is 
approximated  with  the  usual  five-point  scheme.  For  a  cartesian  grid  in  two 
dimensions  the  determinant  of  the  symbol  is 


r4h"2(sin2  |  C,  +  sin2  |  C2> 


d(  ) 


det 


4h"2(sin2  — ■  C,  +  sin2  x  C„) 


2  1 


2  ’2' 


d(t2) 


-  dtC^ 


-  d(C2) 


-  4h”2 ( sin2  ^  C1  +  sin2  j  «2>( |d(C1> |2  +  |d(C2> |2> 


where 


d(C)  -  ih_1sin  C  -  ah”^  iC(2i  sin  j  C)3 


-  2ih-1sin  j  C(cos  J  C  +  4oe1/2  iC 


.  2  1 
sin  j  C) 


Since  d(C)  is  not  zero  for  any  value  of  C,  when  a  is  non-zero,  the 
scheme  is  regular.  Note  that  for  a  equal  to  one-sixth  the  approximations 
(3.2)  and  (3.3)  are  third-order  accurate. 


Since  the  regularized  central  difference  scheme  is  a  variant  of  the 
central  difference  scheme  it  is  easy  to  implement  with  coordinate  maps.  At 
those  boundary  points  where  the  correction  term  would  require  points  beyond 
the  boundary  we  use  the  correction  term  which  interchanges  the  forward  and 
backward  operators.  This  scheme  also  requires  the  use  of  extrapolation  to 


compute  the  pressure  values  on  the  boundary.  It  has  been  found  that  third 
order  extrapolation  gave  quite  good  results,  e.g. 


<3*4)  p0j  “  3p1j  "  3p2 j  +  p3j 

at  the  boundary  x  ■  0  in  two  dimensions. 

A  number  of  first-order  accurate  schemes  for  the  Stokes  and  Navier-Stokes 
equations  have  been  presented  e.g.  Kzivickii  and  Ladyzhenskaya  (1966)  and 
Temam  (1979,  p.  48).  In  this  paper  we  are  concerned  only  with  second-order 
accurate  schemes. 

4.  THE  INTEGRABILITY  CONDITION.  Each  of  the  schemes  for  the  Stokes 
equations  which  have  been  discussed  in  the  previous  section  can  be  written  as 

a)  Vh  +  Vh  -  *h 

(4.1) 

b)  h-k  ■  *h 

with  Dirichlet  boundary  conditions 

km\  on  aah  • 

♦  ♦ 

The  differnce  operators  L  ,  G.  ,  and  Dh  are  approximations  to  the  + 

differential  operators  £n  (1.1T.  £he  discrete  functions  fh,  g h,  and  bh 

are  approximations  to  f,  g  and  b  on  the  mesh  Q  ,  where  hn  is  some  n 

measure  of  the  fineness  of  the  mesh  0.  . 

h 

Now  Jet  us  compare  the  system  (4.1)  with  the  system  (1.1).  First  note 
that  if  is  a  consistent  approximation  to  the  gradient  then  the  discrete 

pressure  p^  is  determined  only  up  to  a  constant.  This  means  that  the  system 
of  linear  equations  (4.1)  does  not  have  full  column  rank.  If  there  are  as 
many  equations  in  (4.1)  as  there  are  unknowns,  and  this  is  the  case  for  each 
scheme  we've  considered,  then  the  system  (4.1)  does  not  have  full  row  rank 
either.  This  implies  that  there  is  a  constraint  which  the  data  must  satisfy 
to  guarantee  a  solution,  in  particular,  the  discrete  integrability  condition 
analogous  to  (1.2)  must  be  satisfied. 

There  are  at  least  two  ways  to  satisfy  the  discrete  integrability 
condition.  The  first  method  would  be  to  analyze  the  matrix  corresponding  to 
(4.1)  and  determine  the  null  space  of  the  adjoint  matrix.  If  the  data  is 
constrained  to  be  orthogonal  to  this  null  space  then  a  solution  will  exist. 
This  approach  is  impractical  for  many  situations  especially  if  coordinate 
changes  have  been  employed  since  then  the  matrices  are  not  easy  to  analyze. 

A  second  approach,  which  will  be  adopted  here,  is  to  replace  (4.1b)  by 

♦  8h 

where  6.  is  a  constant  chosen  to  guarantee  a  solution.  The  value  of  6 


must  be  determined  as  part  of  the  solution.  As  shown  in  the  example:  in 
Section  56  is  at  least  0(h2)  for  the  regularized  central  scheme.  We 
will  refer  to  the  equations  (4.1a(  b* ,  c)  as  (4.1*). 

It  is  interesting  to  note  that  for  the  staggered  mesh  scheme  on  a  uniform 
grid  one  can  easily  satisfy  the  discrete  integrability  condition  since  the 
calculus  of  finite  differences  mimics  the  differential  calculus  very  closely, 
see  e.g.  Kzivickii  and  Ladyzhenskaya  (1966).  Also,  Ghia,  Hankey,  and  Hodge 
(1977)  mention  being  unable  to  obtain  a  solution  to  the  discrete  Navier-Stokes 
equations  for  certain  situations.  We  conjecture  that  this  difficulty  was 
caused  by  the  discrete  integrability  condition  not  being  satisfied. 

There  is  the  possibility  that  the  null  space  of  the  discrete  operator  of 

(4.1)  has  dimension  greater  than  one.  The  regularized  central  scheme  with  the 
third-order  extrapolation  (3.4)  appears  to  have  only  a  one-dimensional  null 
space.  However,  for  a  equal  to  zero  numerical  experiments  indicate  that 
there  are  solutions  which  are  effectively  null  vectors  in  that  they  solve 

(4.1)  with  f  and  g^  smaller  than  the  norm  of  the  solution  by  a  factor 
proportional  to  h  or  h2.  The  dimension  of  the  space  of  nearly  null  vectors 
and  null  vectors  appears  to  be  four  for  the  central  differencing  scheme. 

These  vectors  correspond  to  the  four  zeroes  of  the  determinant  of  the  symbol 
of  the  difference  operator. 

These  nearly  null  vectors  and  null  vectors,  other  than  the  usual  constant 
pressure  null  solution,  make  solving  the  discrete  system  very  difficult.  On 
the  other  hand  the  regular  discrete  systems  can  be  solved  easily  by  the 
iterative  procedure  given  in  Strikwerda  (1982). 

5.  COMPUTATIONAL  RESULTS.  In  this  section  we  present  the  results  of 
testing  the  new  scheme  described  in  Section  3.  In  the  examples  discussed  here 
the  discrete  Stokes  equations  were  solved  using  test  problems  which  illustrate 
various  features  of  the  schemes.  For  each  example  an  exact  analytical 
solution  is  known  and  the  approximate  solutions  were  compared  to  the  exact 
solutions  to  study  the  accuracy  of  the  method.  The  value  of  a,  the 
regularity  parameter,  was  one-sixth  in  all  cases. 

For  the  first  test  problem  the  Stokes  equations  were  solved  on  the  unit 
square  with  a  uniform  grid.  The  exact  soltion  is 

u  =  (2ir)  ^  sin  wx  cos  *y 

(5.1)  v  *  (2x)  ^  cos  wx  sin  *y 

p  *  cos  *x  cos  *y 

♦ 

with  f  *  0  and  g  =  cos  *x  cos  *y.  For  this  example  both  the  accuracy  and 
symmetry  of  the  solution  were  checked.  The  symmetry  was  checked  to  study  the 
effect  of  the  nonsymmetric  regularizing  term  on  the  symmetry  of  the 
solution.  The  symmetry  was  measured  by  computing  the  quantities  sym(u)  and 
sym(p)  given  by 


for  the  third  equation.  The  regularizing  terms  were  added  only  to  the  terms 
corresponding  to  px  in  the  firBt  equation,  py  in  the  second,  and  ux 
and  Vy  in  the  third. 

In  the  third  test  problem  the  Stokes  equations  were  solved  on  a  disc 
using  polar  coordinates  with  uneven  grid  spacing  in  both  the  radial  and 
angular  direction.  The  exact  solution  is 
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r  sin  20 


v  ■  2r  cos  26 


p  *  6r  sin  20 

♦ 

with  £  and  g  being  zero.  The  uneven  grid  was  given  by 

ri  “  *75  pi  +  *25  Pi 


0^  *  ifij  -  .25  sin  <fi 

where  p  and  *f> .  were  evenly  spaced  in  the  interval  [0,1]  and  [0,2v] 
respectively.  This  uneven  spacing  was  chosen  merely  to  show  the  versatility 
o£  the  scheme  and  is  not  intended  to  give  a  better  resolution  of  the  solution. 

For  completeness  we  give  the  Stokes  equations  in  polar  coordinates 
r_1(rur)r  +  r‘2u00  -  r"2u  -  2r"2v0  -  Pr  -  0 
(5.5)  r  ’(rv^  ♦  r  2v00  -  r  2v  +  2r  2u0  -  r  *p0  «  0 

r_1(ru)r  +  r'Vg  =  0  . 


The  difference  formulas  used  in  the  numerical  experiments  we^e  all  second 
order  accurate.  As  an  example  of  the  formulas,  the  term  r  (ru  )  was 
differenced  as 


-  v3>  -  -  v,.,)'  1  'v.  -  i,)  • 


The  results  of  the  numerical  experiments  are  shown  in  the  following 
tables.  Each  table  lists  the  errors  incurred  for  grids  with  N  +  1  points  on 
a  side  for  values  of  N  of  20,  30,  40  and  60.  Tables  1,  II  and  III  list  the 
relative  errors  for  test  problems  1,2  and  3,  respectively,  and  Table  I  also 
shows  £he  symmetry  errors  for  problem  1.  The  relative  errors  are  measured  in 
the  i  -norm  i.e. 


err(u)  -  (£  (u^  -  “(x^y^  )2]f2  /lul2  . 

Also  shown  is  the  value  of  6  which  is  described  in  Section  4.  Table  IV 
displays  the  behavior  of  the  error  as  the  grid  resolution  is  increased.  The 
numbers  shown  are  values  of 


logter^/err^ 

log(N1/N2) 
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TABLE  I 


N 

err(u) 

err(p)  6. 

sym(u) 

sym(p) 

20 

.35 (-3) 

. 1 7 ( —2 )  -.44 (-5) 

.68 (-3) 

. 13 (-2 ) 

30 

.11(-3) 

. 86 ( —3 )  -.89(-6 ) 

.22(-3 ) 

.37 (-3 ) 

40 

.41 (-4) 

.51 (-3)  -.53 (-6) 

.82 (-4) 

. 1 5 ( —3 ) 

60 

. 1 9 ( -4 ) 

•23{-3)  -.50(-7) 

.37 (-4) 

.52 (-4) 

Errors 

values 
-35 (-3) 

for 

of 

test  problem  1  for  grids  with  N  +  1  points  on  a  side  for 
N.  The  numbers  In  parenthesis  are  the  decimal  exponents  1. 
.35  x  IQ-3. 

3 

0  • 

0) 

TABLE  11 

N 

err(u) 

err(p) 

5h 

20 

. 10 (-3 ) 

.2K-2) 

-.24(-3) 

30 

.45 ( —4 ) 

.92 (-3) 

-. 12 (-3 ) 

40 

.25 (-4) 

. 48(-3 ) 

-.74 (—4 ) 

60 

.11(-4) 

.22(-3) 

-•35(-4) 

Errors  for  test  problem  2. 

Table  III 

N 

err(u) 

err(p) 

ih 

20 

. 75  <  —  1 ) 

• 93 ( — 1 ) 

-.33(-2) 

30 

«33(-1 ) 

.34  C — 1 ) 

-.53 (-3) 

40 

.  1 9  <  —  1 ) 

. 18 ( —1 ) 

-. 1 5  <  —3 ) 

60 

•83(-2 ) 

•75(-2 ) 

-.27(-4) 

Errors  for  test  problem  3. 
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TABLE  TV 


N,/N2 

1 

2 

3 

u 

2.8 

2.0 

2.0 

30/20 

P 

1.7 

2.0 

2.5 

& 

h 

4.0 

1.7 

4.5 

u 

3.4 

2.0 

1.9 

40/30 

P 

1.8 

2.3 

2.2 

6h 

1.9 

1.7 

4.4 

u 

3.1 

2.0 

2.0 

40/20 

P 

1.7 

2.1 

2.4 

6 

h 

3.1 

1.7 

4.5 

u 

2.5 

2.0 

2.0 

60/30 

P 

1.9 

2.1 

2.2 

5 

h 

4.2 

1.8 

4.3 

u 

1.9 

2.0 

2.0 

60/40 

p 

2.0 

1.9 

2.2 

5.8 

1.8 

4.2 

Computed  order  of  accuracy  for  u,  p, 


and  <5  for  the  test  problems, 
n 
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where  err^  and  errj  are  the  errors  for  grids  of  N^  +  1  and  Nj  +  1 
points  on  a  side,  respectively.  This  value  should  be  approximately  2.0  for  a 
second-order  scheme.  The  error  reductions  are  shown  for  u,  p  and  6  .  The 
other  velocity  component  had  a  similar  error  behavior  in  all  the  examples. 

All  of  the  solutions  were  computed  by  the  iterative  method  given  in  Strikwerda 
(1982). 

That  some  of  the  errors  were  better  than  second-order  accurate  for  test 
problems  1  and  3  can  be  attributed  to  the  third-order  accurate  difference 
formulas  used  for  the  gradient  and  divergence  terms.  One  might  expect  that 
some  of  the  errors  would  behave  as  third-order  errors  for  some  value  of 
and  n2«  However,  since  the  discrete  Laplacian  is  second-order  accurate, 
for  N  large  enough  the  total  scheme  should  be  second-order  accurate.  It  is 
not  clear  why  <5^  should  behave  as  a  fourth-order  error  as  seen  in  test 
problem  3  and  for  some  values  of  and  N2  in  test  problem  1.  Test 

problem  2  was  no  better  than  second-order  accurate  since  the  gradient  and 
divergence  were  only  second-order  accurate.  The  third-order  differences  were 
only  used  on  those  terms  which  were  necessary  to  achieve  regularity  of  the 
scheme.  The  results  show  conclusively  that  the  scheme  has  overall  second- 
order  accuracy. 

6.  CONCLUSION.  In  this  paper  we  have  examined  several  finite  difference 
methods  for  the  steady  Stokes  and  incompressible  Navier-Stokes  equations  in 
primitive  variables.  We  have  shown  that  the  regularized  centered  difference 
scheme  is  second-order  accurate  and  useful  with  non-rectangular  regions. 
Although  the  numerical  experiments  were  done  using  the  Stokes  equations,  for 
which  exact  solutions  were  available,  we  believe  the  regularized  central 
scheme  is  equally  useful  with  the  incompressible  Navier-Stokes  equations  at 
moderate  Reynolds  number. 
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